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The random-phase approximation (RPA) as an approach for computing the electronic correlation 
energy is reviewed. After a brief account of its basic concept and historical development, the paper 
is devoted to the theoretical formulations of RPA, and its applications to realistic systems. With 
several illustrating applications, we discuss the implications of RPA for computational chemistry 
and materials science. The computational cost of RPA is also addressed which is critical for its 
widespread use in future applications. In addition, current correction schemes going beyond RPA 
and directions of further development will be discussed. 



I. INTRODUCTION 

Computational materials science has developed into an 
indispensable discipline com.plementary to experimental 
materials science. The fundamental aim of computa- 
tional materials science is to derive understanding en- 
tirely from the basic laws of physics, i.e., quantum me- 
chanical first principles, and increasingly also to make 
predictions of new properties or new materials for specific 
tasks. The rapid increase in available computer power to- 
gether with new methodological developments are major 
factors in the growing impact of this field for practical 
applications to real materials. 

Density functional theory (DFT) [l| has shaped the 
realm of first principles materials science like no other 
method today. This success has been facilitated by the 
computational efficiency of the local-density [2| or gener- 
alized gradient approximation [S-Q (LDA and GGA) of 
the exchange-correlation functional that make DFT ap- 
plicable to polyatomic systems containing up to several 
thousand atoms. However, these approximations are sub- 
ject to several well-known deficiencies. In the quest for 
finding an "optimal" electronic structure method, that 
combines accuracy and tractability with transferability 
across different chemical environments and dimension- 
alities (e.g. molecules/clusters, wires/tubes, surfaces, 
solids) many new approaches, improvements and refine- 
ments have been proposed over the years. These have 
been classified by Perdew in his "Jacob's ladder" hierar- 
chyi. 

In this context, the treatment of exchange and cor- 
relation in terms of "exact-exchange plus correlation in 
the random-phase approximation" [7|, |8| offers a promis- 
ing avenue. This is largely due to three attractive fea- 
tures. The exact-exchange energy cancels the spurious 
self-interaction error present in the Hartree energy ex- 
actly (although the RPA correlation itself does contain 
some "self-correlation" and is non-zero for one-electron 
systems). The RPA correlation energy is fully non- 
local and includes long-range van der Waals (vdW) in- 
teractions automatically and seamlessly. Moreover, dy- 
namic electronic screening is taken into account by sum- 
ming up a sequence of "ring" diagrams to infinite order. 



which makes RPA applicable to small-gap or metallic sys- 
tems where finite-order many-body perturbation theories 
breakdown d-dO]. 

The random-phase approximation actually predates 
density-functional theory, but it took until the late 1970s 
to be formulated in the context of DFT [ll| and until 
the early years of this millennium to be applied as a first 
principles electronic structure method [12, [l3 ■ We take 
the renewed and widespread interest of the RPA [la - l4a | 
as motivation for this review article. To illustrate the 
unique development of this powerful physical concept, 
we will put the RPA into its historical context before 
reviewing the basic theory. A summary of recent RPA 
results demonstrates the strength of this approach, but 
also its current limitations. In addition we will discuss 
some of the most recent schemes going beyond RPA, in- 
cluding rcnormalizcd second order perturbation theory 
(r2PT), which is particularly promising in our opinion. 
We will also address the issue of computational efficiency 
which, at present, impedes the widespread use of RPA, 
and indicate directions for further development. 



A. Early history 

During the 1950s, quantum many-body theory under- 
went a major transformation, as concepts and techniques 
originating from quantum electrodynamics (QED) — in 
particular Feynman-Dyson diagrammatic perturbation 
theory — were extended to the study of solids and nu- 
clei. A particularly important contribution at an early 
stage of this development was the RPA, a technique in- 
troduced by Bohm and Pines in a series of papers pub- 
lished in 1951-1953 (HH-Iii]. In recent years, the RPA 
has gained importance well beyond its initial realm of 
application, in computational condensed-matter physics, 
materials science, and quantum chemistry. While the 
RPA is commonly used within its diagrammatic formu- 
lation given by Gell-Mann and Brueckner [ig] , it is nev- 
ertheless instructive to briefly discuss the history of its 
original formulation by Bohm and Pines. Some of the 
cited references are reprinted in [47[ , which also recounts 
the history of the RPA until the early 1960s. Historical 



accounts of the work of Bohm and Pines can be found in 
Refs. [ii-i^ as well as in Refs. ^\^. 

In 1933-1934, Wigner and Seitz published two papers 
on the band structure of metallic sodium |5ll.l52l| in which 
they stressed the importance of the electronic correlation 
energy correction to band-theory calculations of the co- 
hesive energy of metals. Wigner subsequently studied 
the interaction of electrons in a homogeneous electron 
gas (HEG) within a variational approach going beyond 
Hartree-Fock [53 . He initially provided estimates for the 
correlation energy only in the low- and high-density lim- 
its, but later interpolated to intermediate densities |54| . 
For almost two decades, Wigner's estimate remained the 
state of the art in the prototypical many-body problem of 
the HEG. According to a later statement by Herring, "the 
magnitude and role of correlation energy remained inad- 
equately understood in a considerable part of the solid- 
state community for many years." |55l . pp. 71-72] 

Due to the long-range nature of the Coulomb inter- 
action and the resulting divergences, perturbative ap- 
proaches, so successful in other areas, had to be com- 
plemented with approximations that accounted for the 
screening of the charge of an electron by the other elec- 
trons. Before the 1950s, these approximations commonly 
drew on work on classical electrolytes by Debye and 
Hiickel [5a, [531 and on work on heavy atoms by Thomas 
[58| and Fermi [5^,[60l , as well as later extensions J6lll63 |. 

As a reaction to work by Landsberg and Wohlfarth 
[631 164| . Bohm and Pines in 1950 reported to have been 
led "independently to the concept of an effective screened 
Coulomb force as a result of a systematical classical and 
quantum-mechanical investigation of the interaction of 
charges in an electron gas of high density" (65l . p. 103]. 
Their 1951-1953 series of papers [tI. I43l - l45| presents this 
systematical investigation. The RPA was one of several 
physically-motivated approximations in the treatment of 
the HEG which allowed them to separate collective de- 
grees of freedom (plasma oscillations) from single-particle 
degrees of freedom (which today would be called quasi- 
particles or charged excitations) via a suitable canonical 
transformation reminiscent of early work in QED 66 68| . 
A similar theory was developed rather independently for 
nuclei by Bohr and Mottelson (69j . 

In their first paper, illustrating the fundamental idea 
of separating single-particle and collective degrees of free- 
dom, Bohm and Pines introduce RPA as one of four re- 
quirements [43| : 

"(3) We distinguish between two kinds of 
response of the electrons to a wave. One 
of these is in phase with the wave, so that 
the phase difference between the particle re- 
sponse and the wave producing it is indepen- 
dent of the position of the particle. This is 
the response which contributes to the orga- 
nized behavior of the system. The other re- 
sponse has a phase difference with the wave 
producing it which depends on the position of 
the particle. Because of the general random 



location of the particles, this second response 
tends to average out to zero when we consider 
a large number of electrons, and we shall ne- 
glect the contributions arising from this. This 
procedure we call the random phase approxi- 
mation.^^ 



In their second paper [4J], Bohm and Pines develop a 
detailed physical picture for the electronic behavior in 
a HEG due to the presence of Coulomb interactions. 
Only in their third paper, Bohm and Pines treat the 
(Coulomb-) interacting HEG quantum- mechanically ^. 
The RPA enables Bohm and Pines to absorb the long- 
range Coulomb interactions into the collective behavior 
of the system, leaving the single-particle degrees of free- 
dom interacting only via a short-range screened force. 
The RPA amounts to neglecting the interaction between 
the collective and the single-particle degrees of freedom. 
Consequently, the momentum transfers of the Coulomb 
potential in Fourier space can be treated independently. 
The fourth paper [45| applies the new method to the 
electron gas in metals, discussing both validity and con- 
sequences of the RPA, such as the increase in electronic 
effective mass. 

Within condensed-matter theory, the significance of 
the Bohm-Pines approach quickly became apparent: 
Renormalizing the long-range Coulomb interaction into 
an effective screened interaction between new, effective 
single-particle degrees of freedom allowed both to over- 
come the divergences appearing in older theories of in- 
teracting many-body systems and to explain the hith- 
erto puzzling success of the single-particle models of early 
condensed- matter theory (see, e.g., [M])- An early appli- 
cation of the RPA was Lindhard's calculation of the di- 
electric function of the electron gas [701 ■ Alternative ap- 
proaches to and extensions of the Bohm-Pines approach 
were formulated by Tomonaga [7l|, UM : ^^id by Mott [Tj] , 



Frohlich and Pelzer [7^, and Hubbard [75[, 7£|. 

In 1956, Landau's Fermi liquid theory [77| delivered 
the foundation for effective theories describing many- 
body systems in terms of quasiparticles. Brueckner [TSJ 
already in 1955 had introduced a "linked-cluster ex- 
pansion" for the treatment of nuclear matter (sec also 
Ref. [zi]). In 1957, Goldstone (H, using Feynman-like 
diagrams (based on Ref. [8l[), was able to show that 
Brueckner's theory is exact for the ground-state energy 
of an interacting many-fermion system. This put the 
analogy between the QED vacuum and the ground state 
of a many-body system on firm ground. It had been in- 
troduced explicitly by Miyazawa for nuclei [82| and by 
Salam for superconductors [SJ], although the essence of 
the analogy dated back to the early days of quantum field 
theory in the 1930s. 

In late 1956, Gell-Mann and Brueckner employed a 
diagrammatic approach for treating the problem of the 
interacting electron gas. Their famous 1957 paper [4g 
eliminated the spurious divergences appearing in previ- 
ous approaches. Expressing the perturbation scries for 
the correlation energy of the HEG in terms of the Wigner- 



Seitz radius rg, they found that the divergences within 
earher calculations [e.g.,[8Jl were mere artifacts: The log- 
arithmic divergence appearing in the perturbative expan- 
sion of the correlation energy is canceled by similar di- 
vergences in higher-order terms. Summing the diagrams 
(which had a ring structure) to infinite order yielded a ge- 
ometric series, and thus a convergent result. Gell-Mann 
and Brucckner derived an expression for the ground state 
energy of the interacting electron gas in the high-density 
limit. Their work, and Goldstone's paper [80|, are the 
earliest examples of the application of Feynman-type di- 
agrammatic methods in condensed-matter theory. 

Many applications of the new quantum-field theoreti- 
cal methods followed: Gell-Mann calculated the specific 
heat of the high-density HEG [Hj; Hubbard Slla pro- 
vided a description of the collective modes in terms of 
maiiy-body perturbation theory (MBPT); Sawada et al. 
[88l l89j demonstrated that the GcU-Mann-Brucckner ap- 
proach indeed contained the plasma oscillations of Bohm 
and Pines Q , a point around which there had been quite 
some confusion initially [see, e.g.,|90[. In addition, they 
demonstrated that the RPA is exact in the high-density 
limit. In 1958, Nozieres and Pines formulated a many- 
body theory of the dielectric constant and showed the 
equivalence of Gell-Mann and Brueckner's diagrammatic 
approach and the RPA [9l|,|9l|. 



B. RPA in modern times 

Today, the concept of RPA has gone far beyond the 
domain of the HEG, and has gained considerable im- 
portance in computational physics and quantum chem- 
istry. As a key example, RPA can be introduced within 
the framework of DFT [l| via the so-called adiabatic- 
connection fluctuation-dissipation (AGED) theorem [111, 
|93| . |94| . Within this formulation, the unknown exact 
exchange-correlation (XG) energy functional in Kohn- 
Sham [2] DET can be formally constructed by adiabati- 
cally switching on the Goulomb interaction between elec- 
trons, while keeping the electron density fixed at its phys- 
ical value. This is formulated by a coupling-strength in- 
tegration under which the integrand is related to the lin- 
ear density-response function of fictitious systems with 
scaled Goulomb interaction. Thus, an approximation to 
the response function directly translates into an approxi- 
mate DET XG energy functional. RPA in this context is 
known as an orbital-dependent energy functional [95| ob- 
tained by applying the time-dependent Hartree approxi- 
mation to the density response function. 

The versatility of RPA becomes apparent when con- 
sidering alternative formulations. Eor instance, the RPA 
correlation energy may be understood as the shift of the 
zero-point plasmon excitation energies between the non- 
interacting and the fully interacting system, as shown by 
Sawada for the HEG [S^] , and derived in detail by Eurchc 
[96| for general cases (see also Ref. [31[). In quantum 
chemistry, RPA can also be interpreted as an approxi- 



mation to coupled-cluster doubles (GGD) theory where 
only diagrams of "ring" structure are kept [l^ |93| . The 
equivalence of the "plasmon" and "ring-GGD formula- 
tion" of RPA has recently been established by Scuseria 
et al. |15| . These new perspectives not only offer more 
insight into the theory, but also help to devise more ef- 
ficient algorithms to reduce the computational cost, e.g., 
by applying the Gholcsky decomposition to the "ring- 
GGD" equations [l|. 

Following the early work on the HEG, other model 
electron systems were investigated, including the HEG 
surfa ce |98l . |99| . jellium slabs 2J] and jellium spheres 
jlOOJ I . The long-range behavior of RPA for spatially well- 
separated closed-s hell subsystems was examin ed by Sz - 
abo and Ostlund [lOl|, as weh as by Dobson [l0l^U03 |. 
These authors showed that RPA yields the correct 1/R^ 
asymptotic behavior for the subsystem interaction. In 
addition, the long-range dispersion interaction of RPA 
is fully consistent with the monomer polarizability com- 
puted at the same level of theory. This is one of the 
main reasons for the revival of the RPA in recent years, 
because this long-range interaction is absent from LDA, 
GGA, and other popular density-functionals. Other rea- 
sons are the compatibility of the RPA correlation with 
exact exchange (which implies the exact cancellation of 
the self-interaction error present in the Hartree term) and 
the applicability to metallic systems. 

For the HEG it has been demonstr ated that RPA is not 
accurate for short-range correlation |l05l Il06l | , and hence 
for a long time RPA was not considered to be valuable 
for realistic s ystem s. Perdew and coworkers investigated 
this issue [9^, Il07l | , and found that a local/scmi-local cor- 
rection to RPA has little effect on is o- electronic energy 
differences, which suggests that RPA might be accurate 
enough for many practical purposes. The application of 
RPA to realistic systems appeared slightly later, starting 
with the pioneering work of Furche [12| , and Fuchs and 
Gonze [I3 for small molecules. Accurate RPA total en- 
ergies for closed-shell atoms were obtained by Jiang and 
Engel [201 . Several groups investigated molecular proper- 
ties, in particular in the weakly bound r egini e , wit h RPA 
and its variants [ll [13, SI, llfl IH, H, Imjioi, while 
others applied RPA to periodic systems [II [Ull [TTol - 
Ill2l |. Harl and Kresse in particular have performed ex- 
tensive RPA benchmark stud ies for crystalline solids of 
all bonding types [2ll l23 . Ill2| . At the same time, the ap- 
plication of RPA to s urface ad sorption problems has been 
reported |25l - [27l . 133 . Ill3l Ill4l | , with considerable success 
in resolving the "GO adsorption puzzle" . 

Most practical RPA calculations in recent years have 
been performed non-sclf-consistcntly based on a preced- 
ing LDA or GGA reference calculation. In these calcu- 
lations, the Goulomb integrals are usually not antisym- 
metrized in the evaluation of the RPA correlation energy, 
a practice sometimes called direct RPA in the quantum- 
chemical literature. In this paper we will denote this com- 
mon procedure "standard RPA" to distinguish it from 
more sophisticated procedures. While a critical assess- 



mcnt of RPA is emerging and a wide variety of applica- 
tions are pursued, certain shortcomings of standard RPA 
have been noted. The most prominent is it s sy stematic 
underestimation of binding energies I13 . l2a 1341 1 , and the 



failure to describe stretched radicals [l8|, Ill5l llla | . Over 



the years several attempts have been made to improve 
upon RPA. The earliest is RPA-I-, where, as mentioned 
above, a local/semi-local correlation correction based on 
LDA or GGA is add ed to the standard RPA correla- 
tion energy [9^, Il07l |. Based on the observation that 
in molecules the correlation hole is not sufficiently ac- 
curate at medium range in RPA, this has recently been 
extended to a non-local correction sche me |36l . |37| . Sim- 
ilarly, range-separated frameworks J117l | have been tried, 
in which onl y th e l ong-range part of RPA is explicitly 
included |16l. IitI. l28l . Ill8l4l22l |. whereas short/mid-range 
correlation is treated differently. Omitting the short- 
range part in RPA is also numerically beneficial whereby 
the slow convergence with respect to the number of basis 
functions can be circumvented. Due to this additional 
appealing fact, range-separated RPA is now an active 
research domain despite the empirical parameters that 
govern the range separation. Another route to improve 
RPA in the framework of ACFD is to add an /xc kernel 
to the response fu nction an d to find suitable approxima- 
tions for it [ll|3i,[l2l[l2l. Last but not least, the COD 
perspective offers a different correction in form of the 
second -orde r screened exchange (SOSEX) contribution 
la . [971 [l25| ■ whereas the MBPT perspective inspired sin- 
gle excitation (SE) corrections [33|. SOSEX and SE are 
disti nct many-body corrections and can also be combined 
[l26J | . These corrections have a clear diagrammatic repre- 
sentation and alleviate the above-mention ed un derbind- 
ing problem of standard RPA considerably |l26l | . Yet an- 
other proposal to improve RPA by incorporating higher- 
order exchange effect s in various ways has also been dis- 
cussed recently 121|. However, at this point in time, a 
consensus regarding the "optimal" correction that com- 
bines both efficiency and accuracy has not been reached. 

Although the majority of practical RPA calculations 
arc performed as post-processing of a preceding DFT 
calculation, self-consistent RPA calculations have also 
been performed within the optimized effective potential 
(OEP) framework. OEP is a procedure to find the opti- 
mal local multiplicative potential that minimizes orbital- 
dependent energy functionals. The first RPA-OEP cal- 
culations actually date back more than 20 years, but 
were not recognized as such. Godby, Schliiter and Sham 
solved the Sham- Schliiter equation for the GW self- 
energy [l27l.[T2i|. which is equivalent to the RPA-OEP 
equation, for the self-consistent RPA KS potential of bulk 
silicon and other semiconductors, but did not calculate 
RPA ground-state energies. Similar calculations for other 
bulk materials followed later by Kotani |l29| and Griining 
et al. |l30l . Il3l| . Hellgren and von Earth 132| and then 
later Verma and Bartlett [ij have looked at closed-shell 
atoms and observed that the OEP-RPA KS potential 
there reproduces the exact asymptotic behavior in the 



valence region, although its behavior near the nucleus 
is not very accurate. Extensions to diatomic molecules 
have also appeared recently [4l|, ^^. Our own work 
on the SE correction to RPA indicates that the input- 
orbital dependence in RPA post-processing calculations 
is a significant issue. Some form of self-consistency would 
therefore be desirable. However, due to the considerable 
numerical effort associated with OEP-RPA calculations, 
practical RPA calculations will probably remain of the 
post-processing type in the near future. 



Despite RPA's appealing features its widespread use in 
chemistry and materials science is impeded by its com- 
putational cost, which is considerable compared to con- 
ventional (scmi)local DFT functionals. Furche's original 
implementation based on a molecular particle-hole basis 
scales as 0{N^) [ij]. This can be reduced to 0{N^) using 
the plasmon-pole formulation of RPA [9al, or to 0{N^) 
[201 when the resolution-of-identity (RI) technique is em- 
ployed. Scuseria et al. [15[ pointed out even slightly ear- 
lier that a 0{N^) scaling can be achieved by combining 
the "ring-CCD" RPA formulation and the Cholesky de- 
composition techn ique . Our own RPA implementation 
[l33| in FHI-aims |l34| . whi ch has been used in produc- 
tion calculations |26l. l34lll26| before, is based on localized 
numeric atom-centered orbitals and the RI technique, 
and hence naturally scales as 0{N^). Plane- wave based 
implementations 13l . lll2J also automatically have ©(Af'*) 
scaling. However, in standard implementations the con- 
vergence with respect to unoccupied states is slow. Pro- 
posals to e l imin ate the dependence on the unoccupied 
states |l35l . Il36l | in the context of plane wave bases, by 
obtaining the respo nse fu nction from density-functional 
perturbation theory [137| . have not been explored so far 
for local-orbital based implementations. In local orbital 
based approaches the 0(7V^) scaling can certainly be re- 
duced by exploiting matrix spars ity, as demonstrated re- 
cently in the context of GW [l38j | or secon d-order M0ller- 
Plessct perturbation theory (MP2) |l39l |. Also approxi- 
mations to RPA [3l| or effective screening models |l4Cll | 
might significantly improve the scaling and the compu- 
tational efficiency. Recently, RPA has been cast into 
the continuum mechanics formulation of DFT [l4l| with 
cons iderable success in terms of computational efficiency 
[142| . In general, there is still room for improvement, 
which, together with the rapid increase in computer 
power makes us confident that RPA-type approaches will 
become a powerful technique in computational chemistry 
and materials science in the future. It would thus be de- 
sirable, if the material science community would start to 
build up benchmark sets for materials scie nce akin to th e 
ones in quantum chemistry (e.g. G2 [l43J | or S22 MJ]). 
These should include prototypical bulk crystals, surfaces, 
and surface adsorbates and would aid the development of 
RPA-based approaches. 



II. THEORY AND CONCEPTS 



nian Hi (A) in Eq. ([T]) becomes 



RPA can be formulated within different theoretical 
frameworks. One particularly convenient approach to 
derive RPA is the so-called "adiabatic connection (AC)", 
which is a powerful mathematical technique to obtain the 
ground-state total energy of an interacting many-particle 
system. Starting with the AC approach, the interacting 
ground-state energy can be retrieved either by coupling 
to the fluctuation-dissipation theorem in the DFT con- 
text, or by invoking the Green-function based MBPT. 
RPA can be derived within both frameworks. In addi- 
tion, RPA is also intimately linked to the coupled-cluster 
(CC) theory. In this section, we will present the theoret- 
ical aspects of RPA from several different perspectives. 



A. Adiabatic connection 

The ground-state total energy of an interacting many- 
body Hamiltonian can formally be obtained via the AC 
technique, in which a continuous set of coupling-strength 
(A) dependent Hamiltonians is introduced 



N . N 
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= E T^ + lE[-r'(r.)-^''=''*(r,)--"^W]. 

(4) 

In the AC construction of the total energy, we introduce 
the ground-state wave function |$a) for the A-dependent 
system such that 



H{X)\^x)=E{X)\^x)^ 



(5) 



Adopting the normalization condition (^a|*a) = 1, the 
interacting ground-state total energy can then be ob- 
tained with the aid of the Hellmann-Feynman theorem, 
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H{X) = Ho + XHi{X), 



(1) 



where 



that "connect" a reference Hamiltonian Hq = H{X = 0) 
with the target many-body Hamiltonian H = H{X = 1). 
For the electronic systems considered here, H{X) has the 
following form: 
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where N is the number of electrons, w^'^* is a A-dependent 
external potential with v"^-^ (r) = v^^^ (r) being the phys- 
ical external potential of the fully-interacting system. 
Note that in general v'^'^ can be non-local in space for 
A 7^ 1. Hartrce atomic units h = e = rUe = 1 arc used 
here and in the following. The reference Hamiltonian Hq, 
given by Eq. ([2]) for A = 0, is of the mean-field (MF) type, 
i.e., a simple summation over single-particle Hamiltoni- 
ans: 
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In Eq. ©, 



,MF 



is a certain (yet-to-be-specified) mean- 



field potential arising from the electron-electron interac- 
tion. It can be the Hartree-Fock (HF) potential v^^ or 
the Hartrce plus exchange-correlation potential v^^"^ in 
DFT. Given Eq. © and ^, the perturbative Hamilto- 



£;o = S(n) = (*o|i^o|*o) 



(7) 



is the zeroth-order energy. We note that the choice of 
the adiabatic-connection path in Eq. ([S]) is not unique. 
In DFT, the path is chosen such that the electron density 
is kept fixed at its physical value along the way. This im- 
plies a non-trivial (not explicitly known) A-dependence 
of Hi{X). In MBPT, one often chooses a linear connec- 
tion path — Hi{X) — Hi (and hence dHi{X)/dX = 0). 
In this case, a Taylor expansion of |^a) in terms of A in 
Eq. ([6]) leads to standard Rayleigh-Schrodinger pertur- 
bation theory (RSPT) [mSJ . 



B. RPA derived from ACFD 

Here wc briefly introduce the concept of RPA in the 
context of DFT, which serves as the foundation for most 
practical RPA calculations in recent years. In Kohn- 
Sham (KS) DFT, the ground-state total energy for an 
interacting A^-electron system is an (implicit) functional 
of the electron density n(r) and can be conveniently split 
into four terms: 

E[n{r)] = Ts[M^)]+E,^t[n{r)]+EH[n{r)]+E^,[ij,{r)] . 

(8) 
Tg is the kinetic energy of the KS independent-particle 
system, E^^t the energy due to external potentials. 
En the classic Hartree energy, and E^c the exchange- 
correlation energy. In the KS framework, the electron 
density is obtained from the single-particle KS orbitals 



■0i(r) via n(r) = X^i"^'^ l''/'i('')l^- Among the four terms in 
Eq. ([8)), only Ecxt[n{r)] and Eii[n{r)] are explicit func- 
tionals of n(r). Tg is treated exactly in KS-DFT in terms 
of the single-particle orbitals ipi (r) which themselves are 
functionals of n{r). 

All the many-body complexity is contained in the un- 
known XC energy term, which is approximated as an 
explicit functional of n(r) (and its local gradients) in con- 
ventional functionals (LDA and GGAs), and as a func- 
tional of the ^i (r) 's in more advanced functionals (hybrid 
density functionals, RPA, etc.). Different existing ap- 
proximations to -Exc can be classified into a hierarchical 
scheme known as "Jacob's ladder" [6| in DFT. However, 
what if one would like to improve the accuracy of E^c in a 
more systematic way? For this purpose it is illuminating 
to start with the formally exact way of constructing E-^c 
using the AC technique discussed above. As alluded to 
before, in KS-DFT the AC path is chosen such that the 
electron density is kept fixed. Equation ^ for the exact 
ground-state total-energy E = E{\ = 1) then reduces to 

/•I ^ d 

= Eo+- / dX drdr'x 

^^' Ir-r'l ' ^' 



where 



+ I drnir) [vTMr) - vTUr) 



N 

n(r) = \J 5{r — r^) 



(9) 



(10) 



is the electron-density operator, and n(r) = 
(*A|?i(r)|*A) for any < A < 1. 

For the KS reference state \'i>o) (given by the Slater 
determinant of the occupied single-particle KS orbitals 
{-0^(1")}) we obtain 



N r 



Eo = (*o|^ 



1 9 



«A=o(r.) 



l*o) 



= TAAir)]+Jdrn{r)vTMr), 
and thus 
E =T, [7/;,(r)] + f drn{r)vTMr)+ 

Mr) [n{r') - 5{r - r')] 



(11) 



ll'dxjjdrdr'i^. 



*a). 
(12) 



Equating ([8]) and ((T2)) . and noticing 

E^Hr)] = i fdrdr'!p^ (13) 



r — r' 



i?ext[n(r)] = y'rfrn(r)«rJi(r), 



(14) 



one obtains the formally exact expression for the XC en- 
ergy 



E. 



lldX 



^^^^,<(r,r>(r) 

r — r' 



(15) 



Here 



is the formal expression for the so-called XC hole, with 
6n(r) = n{r) — n(r) being the fluctuation of the density 
operator n(r) around its expectation value n(r). Equa- 
tion (|16|1 shows that the XC hole is related to the density- 
density correlation function. In physical terms, it de- 
scribes how the presence of an electron at point r depletes 
the density of all other electrons at another point r'. 

In the second step, the density-density correlations 
(fluctuations) in Eq. p6p are linked to the response 
properties (dissipation) of the system through the zero- 
temperature /Zuctuation-dJssJpaiJOTi theorem (FDT). The 
FDT is a powerful technique in statistical physics. It 
states that the response of a system at thermodynamic 
equilibrium to a small external perturbation is the same 
as its response to the spontaneous internal fluctuations 
in the absence of the perturbation [l46| . The FDT 
is manifested in many physical phenomena and applies 
to both thermo and quantum-mechanical fiuctuations. 
The dielectric formu latio n of the many-body problem by 
Nozieres and Pines jl47| is a key example of the latter. 
In this context, the zero-temperature FDT leads to |148| 

(*A|'5n(r)(5n(r')|*A) = / dujlmx\r , r' , uj) , (17) 

T^ Jo 

where x^(r, r',aj) is the linear density response func- 
tion of the (A-scaled) system. Using Eqs. (|15m7|) and 
z;(r, r') = l/|r — r'|, we arrive at the renowned ACFD 
expression for the XC energy in DFT 



E'xc = - / dX // drdr'v{r,r') x 

/ (iajlmx'^(r,r',a;) — (5(r — r')7i(r) 

= — f dX ffdrdr'v{r,r') x 

dujx^iy^, I"', Jw) — 6{r — r')n(r) 



(18) 



The reason that the above frequency integration can be 
performed along the imaginary axis originates from the 
analytical structure of x'*'(r, r', w) and the fact that it be- 
comes real on the imaginary axis. The ACFD expression 
in Eq. (|18p transforms the problem of computing the XC 
energy to one of computing the response functions of a 
series of fictitious systems along the AC path, which in 
practice have to be approximated as well. 

In this context the random-phase approximation is a 
particularly simple approximation of the response func- 
tion: 

Xrpa(i", r', i^) = X°(r, r', iuj) + 

/ dridr2X°(r,ri,iw)Aw(ri - r2)xRPA(i"2,r',a;). 

(19) 

X° (r, ri , iiS) is the independent-particle response function 
of the KS reference system at A = and is known ex- 
plicitly in terms of the single-particle KS orbitals ipi{r), 
orbital energies e^ and occupation factors fi 

0, , . , ^(/.-/,M(r)^,(r)^*(r')V'.(r') 



•y 



e,; — £.; — lOJ 



(20) 

From equations (fTHj) and ([T9|) . the XC energy in RPA can 
be separated into an exact exchange (EX) and the RPA 
correlation term, 



Et 



RPA _ 77EX I 77RPA 



E. 



E 



(21) 



where 



E^^ = - ff drdr'v{r,r') x 

rfcjx''(r, r', i^) — 5{r — r')n{r) 



C. RPA derived from MBPT 

An alternative to ACFD is to compute the interacting 
ground-state energy by performing an order-by-order ex- 
pansion of Eq. (I6|). To this end, it is common to choose 



a linear AC path, i.e., in Eq. (|4]) v™* 



' + {l-X)v 



MF 



such that 
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Now equation ([6]) 


reduces to 






E^ 


Eo+ f 


1 
rfA(*A 


Hil^x) 



,MF 



(25) 



(26) 



A Taylor expansion of |^a) and a subsequent A inte- 
gration lead to an order-by-order expansion of the in- 
teracting ground-state total energy, e.g., the first-order 
correction to Eq is given by 



E 



(1) - 



dXi'^olHil^o) 


l-ol-ffil^o) 



E 



H 



77EX _ rpMF 



(27) 



where E^ and E^^ are the classic Hartrce and exact ex- 
change energy defined in Eq. p^ and (j^ . respectively. 



£:Mi' = (^olt^iii-l^o) is the "double-counting" term due 
to the MF potential v^^ , which is already included in 
H^. The sum of Eq and the first-order term i?^^' yields 
the Hartree-Fock energy, and all higher-order contribu- 
tions constitute the so-called correlation energy. 

The higher-order terms can be evaluated using the di- 
agrammatic technique developed by Goldstone [80|- For 
instance, the second-order energy in RSPT is given by 



^ /./, II drrfrV*(r)^,(r)«(r, r')r,{r')Mr') 



^^^'=E 



\{<Po\Hl\<Pn) 



(22) 



and 

^RPA 



n>0 
occ unocc 

EE 






|($o|i?i|* 



a\|2 



occ unocc 



En-E. 



(0) 



|($o|gi|<l>°^)p 



ij ab 



En-E 



— / / drdr'v(r,r') x 



(0) 

ij,ab 

(28) 



dio 



f^-^XRPA(i", r',iuj) - xo(r, r', «w) 



2^7o 



dujTi [ln(l - x"(iw)w) + x°(«w)w] 



(23) 



For brevity the following convention 

Ti-[AB] = // drdr'A{r,r')B{r',r) (24) 

has been used in Eq. ([^S)) . 



where |$o) = l^o) is the ground state of the reference 
Hamiltonian Hq, and \^n) for n > correspond to its 

excited states with energy En = ($„|-ffo|$n)- \^n) 
can be classified into singly-excited configurations |$i.a), 
doubly-excited configurations |$ij,a&), etc.. The summa- 
tion in Eq. p8|) terminates at the level of double exci- 
tations. This is because Hi only contains one- and two- 
particle operators, and hence does not couple the ground 
state |$o) to triple and higher-order excitations. We will 
examine the single-excitation contribution in Eq. p8[) in 
detail in Section IIVCI Here it suffices to say that this 
term is zero for the HE reference and therefore is not 
included in MP2. The double-excitation contribution 





E^ 



FIG. 1: Goldstone diagrams for the MP2 correlation energy. 
The two graphs describe respectively the second-order direct 
process, and the second-order exchange process. The upgoing 
solid line represents a particle associated with an unoccupied 
orbital energy ta, the downgoing solid line represents a hole 
associated with an occupied orbital energy ei, and the dashed 
line denotes the bare Coulomb interaction. 



can be further split into two terms, corresponding to the 
second-order direct and exchange energy in MP2, whose 
representation in terms of Goldstone diagrams is depicted 
in Fig. [1] The rules to evaluate Goldstone diagram s can 
be found in the classic book by Szabo and Ostlund |l45j . 
As a side remark, the application of MP2 used to be re- 
stricted to finite systems, and its extension to infinite 
periodic systems has been a big challenge because of its 
0{N^) canonical scaling. In recent years, however, sev- 
eral authors demonstrated that it is in principle feasi- 
ble toapply MP2 to one- or two-dimensional systems 
149l - ll5l| . With t he m ore recent implementations in the 
CRYSCOR code [HI as weU as in the VASP code ^, 
the application of MP2 to three-dimensional crystalline 
solids has become realistic. 

The Goldstone approach is convenient for the lowest 
few orders, but becomes cumbersome or impossible for 
arbitrarily high orders, the evaluation of which is essen- 
tial when an order-by-order perturbation breaks down 
and a "seZeciwe summation to infinite order'''' procedure 
has to be invoked. In this case, it is much more conve- 
nient to express the total energy in terms of the Green 
function a nd th e self-energy, as done, e.g., by Luttinger 
and Ward [l53| . Using the Green function langu age, t he 
ground-state total energy can be expressed as [3, Il53| , 



1 r^ d\ / ^ r°° 



(29) 



1 f d\ f ^ r°° 



(30) 



where G° and G(A) arc single-particle Green functions 
corresponding to the non-interacting Hamiltonian Hq 
and the scaled interacting Hamiltonian H{X), respec- 
tively. I]*(A) and S(A) are the proper (irreducible) and 
improper (reducible) self-energies of the interacting sys- 
tem with interaction strength A. [Proper self-energy dia- 
grams are those which cannot be split into two by cutting 
a single Green function line.] Note that in Eq. (|29|) and 
(|5n)) , the trace convention of Eq. ([M)) is implied. 



The above quantities satisfy the following relationship 

G{iuj,X) = G°{iuj) + G°{iuj)J:{iuj,X)G°{iuj) 

= G^iiijj) + G°(ia;)E*(ia;, A)G(ia;, A). (31) 



From Eq. ^Q the equivalence of Eq. dH]) and Eq. ((50)) 
is obvious. In Eq. (j29p . a perturbation expansion of 
the A-dependent self-energy Yi^iiw) naturally translates 
into a perturbation theory of the ground-state energy. 
In particular, the linear term of I]^(ia;) yields the first- 
order correction to the ground-state energy, i.e., E^^' 
in Eq. (P7)) . All higher-order [n > 2) contributions of 
I]^(ia;), here denoted Ec, define the so-called correla- 
tion energy. In general the correlation energy cannot be 
treated exactly. A popular approximation to Ec is the 
GW approach, which corresponds to a selective summa- 
tion of self-energy diagrams with ring structure to infi- 
nite order, as illustrated in Fig. [2l^a). Multiplying Go to 
the GW self-energy Y,'^^ {ibj) as done in Eq. (^9]) and 
performing the A integration, one obtains the RPA cor- 
relation energy 

(32) 
This illustrates the close connection between RPA and 
the GW approach. A diagrammatic representation of 
^RPA jg s]^Q.^^j^ jjj Fig. [5](b). We emphasize that the dia- 
grams in Fig. [2l^a) and[2ljb) are Feynman diagrams, i.e., 
the arrowed lines should really be interpreted as prop- 
agators, or Green functions. A similar representation 
of^RPA ^^^ |-jg drawn in terms of Goldstone diagrams 
[l45j . as shown in Fig[2jc). However, caution should be 
applied, because the rules for evaluating these diagrams 
arc different (see e.g., Ref. J9l ll45l |). and the prefactors in 
Fig. EJb) are not present in the corresponding Goldstone 
diagrams. The leading term in RPA corresponds to the 
second-order direct term in MP2. 

We note that starting from Eq. (P5)l this procedure 
naturally gives the perturbative RPA correlation energy 
based on any convenient non-interacting reference Hamil- 
tonian i/o, such as Hartree-Fock or local/semi- local KS- 
DFT theory. If one instead starts with Eq. (|30|l and 
applies the GW approximation therein, G(A, iuj) and 
S*(A, iuj) become the self-consistent GW Green function 
and self-energy. As a result the improper self-energy di- 
agrams in Eq. ([29]) . which are neglected in the perturba- 
tive GW approach (known as G^W^ in the literature), 
are introduced and the total energy differs from that of 
the RPA. An in-depth discussion of self-consi stent GW 
and its implications can be found in [154l4l57 |. 



D. Link to coupled cluster theory 

In recent years, RPA has also attracted considerable 
attention in the quantum chemistry community. One key 
reason for this is its intimate relationship with coupled 



(a) 



E?^(A) = y 



(c) 



E^"^ = 



<:>:.. •6^"'^ 



(b) 



^RPA ^ _ 1 , 



+ 









where c^ and c are single-particle creation and annihila- 
tion operators and t°, tf'^, . . . are the so-called CC singles, 
doubles, . . . amplitudes yet to be determined. As before, 
i,j, ... refer to occupied single-particle states, whereas 
a, 6, • • • refer to unoccupied (virtual) ones. Acting with 
Tn on the non-interacting reference state |$o) generates 

yabc- 

ijk- 



ri-order excited configuration denoted 1$"^^.'.'.'}: 



r„|$o) 



i>j>k--- .a>b>c--- 



ahc-- 
ijk-- 



\^i]k-- 



■)• 



(37) 



The next question is how to determine the expansion 
coefhcients i°j'^.'.? The CC many-body wave function 
in Eq. (|33l) has to satisfy the many-body Schrodingcr 
equation. 



iJe^|$) =Ee^\<^). 



(38) 



or 



FIG. 2: Feynman diagrams for the GW self-energy (a), Feyn- 
man diagrams for the RPA correlation energy (b), and Gold- 
stone diagrams for the RPA correlation energy (c). Solid lines 
in (a) and (b) (with thick arrows) represent fermion propa- 
gators G, and those in (c) (with thin arrows) denote parti- 
cle (upgoing line) or hole states (down-going line) without 
frequency dependence. Dashed lines correspond to the bare 
Coulomb interaction v in all graphs. 



cluster (CC) theory, which has been very successful for 
accurately describing both covalent and non-covalent in- 
teractions in molecular systems. To understand this re- 
lationship, we will give a very brief account of the CC 
theory here. More details can for inst ance be found in 
a review paper by Bartlett and Musial jl58| . The essen- 
tial concept of CC builds on the exponential ansatz for 
the many-body wave function VE* for correlated electronic 
systems 



I*) 



1$) 



(33) 



$) is a non-interacting reference state, usually chosen to 
be the HF Slater determinant, and T is a summation of 
excitation operators of different order. 



T = Ti + T2 + T3 + • ■ • + r„ 



(34) 



with Ti, T2, T3, ■ • • being the single, double, and triple 
excitation operators, etc. These operators can be most 
conveniently expressed using the language of second- 
quantization, namely. 



(35) 



Ti = 




T2 = 


ij,ab 


T„ = 


1 \ ^ fabc- 

^ ' ijk--- ,abc--- 



. »+ »t ~t 



'bCc' ■ ■ Ck Cj Ci , 



(36) 



e"^iJe^|$) = £^1$). 



(39) 



By projecting Eq. ((M)) onto the excited configurations 
^iik--- ) ' which have zero overlap with the non- interacting 



ground-state configuration |$), one obtains a set of cou- 

^abc- 

^ijk- 



pled non-linear equations for the CC amplitudes t^l^^" 



/^abc- 



|e-^ife^|$) =0. 



(40) 
self- 



These can be determined by solving Eqs. PO)) 
consistently. 

In analogy to the Goldstone diagrams, equation ([40|) 
can be represent ed p ictorially using diagrams, as illus- 
trated by Cizck [1591 in 1966. In practice, the expan- 
sion of the T operator has to be truncated. One popular 
choice is the C C do ubles (CCD) approximation, or T2 
approximation |159| . that retains only the double exci- 
tation term in Eq. p4p . The graphical representation of 
CCD contains a rich variety of diagrams including ring 
diagrams, ladder diagrams, the mixture of the two, etc. 
If one restricts the choice to the pure rii ig di agrams, as 
practiced in early work on the HEG |97l. Il60l |. the CCD 
equation is reduced to the following simplified form |la | , 



B + AT + TA + TBT = 0. 



(41) 



A, B, T are all matrices of rank A'occ • A^vir with A'occ and 
Avir being the number of occupied and unoccupied single- 
particle states, respectively. Specifically we have Aiajb = 
(ej - ea)S^jSab - {ib\aj), B^ajb = {ij\ab), and T^ajb = tfj, 
where the Dirac notation for the two-electron Coulomb 
repulsion integrals 



{pq\rs) 



dvdv' 



>;(r)V.(r)7^*(r')^s(r') 



(42) 



has been adopted. 

Equati on ([¥ 1]) is mathematically known as the Riccati 
equation |l6l| . Solving this equation yields the ring-CCD 
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amplitudes T'"*-^*-'^, with which the RPA correlation en- 
ergy can be written as 



rRPA ^ ixr fST'CCD) 



1 






rCCD 



n / , -^iaJb-L jb,ia 
ij,ab 



(43) 



The CCD formulation of RPA as given by Eq. (14Tj) and 
(P5|l was shown by Scuscria et al. [l^ to be analyti- 
cally equivalent to the plasmonic formulation of RPA. 
The latter has recently been discussed in detail by Furche 
[4Cl l96j , and hence will not be presented in this review. 
Technically, the solution of the Riccati equation (|4T|) is 
not unique, due to the non-linear nature of the equation. 
One therefore has to make a judicious choice for the ring- 
CCD amplitudes in practical RPA calculations |ll6| . 



as a matrix equation of rank A^aux, which is numerically 
very cheap to evaluate. The dominating step then be- 
comes the build of the matrix form of x" which scales 
as O jN'^ ) . We refer the readers to Appendix [X] and 
Ref. |l33l | for further details. 

Plane- wave based implementations [1^ Ill2| automat- 
ically have 0{N'^) scaling. In a sense the plane- wave 
based RPA implementation is very similar in spirit to the 
local-orbital-based RI-RPA implementation. In the for- 
mer case the plane waves themselves serve as the above- 
mentioned ABFs. 



III. ALGORITHMS AND IMPLEMENTATIONS 



B. Speed-up of RPA with iterative methods 



A. RPA implementations and scaling 

In this section we will briefly review different imple- 
mentations of the RPA approach, since scaling and effi- 
ciency are particularly important for a computationally 
expensive approach like the RPA. Also, for historical rea- 
sons, the theoretical formulation of RPA is often linked 
closely to a certain implementation. Similar to conven- 
tional DPT functionals, implementations of RPA can be 
based on local orbitals (LO), or on plane waves, or on 
(linearized) augmented plane waves (LAPW). LO im- 
plementations have been re porte d for the development 
version of Gaussia n lla . llSl . Ill8l | , t he d evelopment ver- 
sion of Molpro [13, Il23l |. FHI-aims [l3|| and Turbomole 
\VA |3C| . Plane- wave based i mple mentations can be found 
in ABINIT [H, VASP |2l|,[il2|, and Quantum-Espresso 
23lll35 |. An early implementation by Miyake et al. 



The RPA correlation energy in (|23|) can also be rewrit- 
ten as follows. 



was based on LAPW. 

Furche's original implementation uses a molecular 
particle- hole basis and scales as 0{N^) [12[, where N is 
the number of atoms in the system (unit cell) . This can 
be reduced to 0{N^) using t he p lasmon-pole formulation 
of RPA m, or to 0{N^) [sj when the rcsolution-of- 
identity (RI) techn ique is empl oyed. Our own RPA im- 



plementation jl33l | in FHI-aims jl34l | is described in Ap- 
pendix |3 It is based on localized numeric atom-centered 
orbitals and the RI technique, and hence naturally scales 
as 0{N'^). 

The key in the RI-RPA implementation is to expand 
the occupied- virtual orbital pair products 0,*(r)^j(r) ap- 
pearing in Eq. (I20p in terms of a set of auxiliary ba- 
sis functions (ABFs) {P^(r)}. In this way, one can re- 
duce the rank of the matrix representation of x° from 

A^occ * A^vir to A^aux with A^aux < A^occ * A^vir- Here A^aux, 

A^occj and A^vir denote the number of ABFs, and the num- 
bers of occupied and unoccupied [virtual) single-particle 
orbitals respectively. With both x° and the coulomb ker- 
nel V represented in terms of the ABFs, the RPA correla- 
tion energy expression in Eq. (|23|) can be rc-interpretcd 



Wa 



E^^ = -iz duY: [In {s^,i^u.)) + 1 - s°(^c.)] , 



2n 

(44) 
where e^{iuj) is the /ith eigenvalue of the dielectric func- 
tion e{iu!) = 1 — x'^{iuj)v represented in the ABFs. All 
eigenvalues are larger than or equal to 1. From (|44p it 
is clear that eigenvalues which are close to 1 have a van- 
ishing contribution to the correlatio n ene rgy. For a set 
of different materials, Wilson et al. |l36l | observed that 
only a small fraction of the eigenvalues differs signifi- 
cantly from 1, which suggests that the full spectrum of 
e{iLo) is not required for accurate RPA correlation en- 
ergies. This opens up the possibility of computing the 
RPA correlation energy by obtaining the "most signifi- 
cant" eigenvalues of e(iw) (or equivalently x'^(iw)z;) from 
an iterative diagonalization procedure, instead of con- 
structing and diagonalizing the full e(iw) or x''(*w)i' ma- 
trices. In practice this can be conveniently done by re- 
sorting to the linear response techniq ue o f density func- 
tional perturbation theory (DFPT) jl37J and has been 
proposed and implemented in the two independent works 
of Galli and coworkers 123, Il36i . and of Nguyen and 



de Gironcoli |135l | within a pseudopotential plane-wave 
framework. In these (plane- wave based) implementations 
the computational cost is reduced from N^^ NoccNvh- to 
A^pw-i/jA^^cc^cig! where A'pw-x and A^pw-i/j are the numbers 
of plane waves to expand the response function x'^ and 
the single-particle orbitals ip respectively, and A^cig is the 
number of dominant eigenvalues. In this way, although 
the formal scaling is still 0{N'^), one achieves a la rge re- 
duction of the prefactor, said to be 100-1000 [l35| . This 
procedure is in principle applicable to RI-RPA implemen- 
tation in local-orbital basis sets as well, but has, to the 
best of our knowledge, not been reported so far. 
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IV. COMPUTIONAL SCHEMES BEYOND RPA 

In this Section we will give a brief account of the ma- 
jor activities for improving the standard RPA, aiming at 
better accuracy. 



E^ 




+ 



S-^ 



+ 



A. Semi-local and non-local corrections to RPA 

It is generally accepted that long-range interactions are 
well described within R PA, w hereas short-range correla- 
tions are not adequate |l06l |. This deficiency manifests 
itself most clearly in the pair-correlation function of the 
HEG, which spuriously becomes negati ve when t he sepa- 
ration between two electrons gets small Il05l.ll06|. B ased 
on this observation, Perdew and coworkers [93, Il07l | pro- 
posed a semi-local correction to RPA, termed as RPA-I- 



E. 



RPA+ _ j^HPA 



j^GGA _ ^GGA-RPA 



(45) 



where E^'^^ is the GGA correlation energy, and 
]^GGA-BPA j-epresents the random-phase approximation 
within GGA. Thus the difference between E^'^^ and 
E^'^^~^^^ gives a semi-local correction to RPA for in- 
homogeneous systems. As mentioned before in the intro- 
duction, the RPA+ scheme, although conceptually ap- 
pealing, and good for total energies [20|, does not sig- 
nificantly improve the description of energy differences, 
in particular the atomization energies of small molecules 
[13 |. This failure has been attributed to the inaccuracy 
of RPA in describing the multi-center non-locality of the 
correlation hole, which cannot be corrected by semi-local 
corrections of the RPA+ type [3a, l37| ■ A fully non-local 
correction (nlc) to RPA has recently been proposed by 
Ruzsinszky, Perdew, and Gsonka [33|. It takes the fol- 
lowing form 



E. 



nlc 



drn(r) [e^^^{r) 



gGGA-RPA , J. 



where e^'^^(r) and e 



GGA-RPA, 



)][l-aF(/(r)]] 

(46) 
r) are the GGA energy 
density per electron and its approximate value within 
RPA, respectively, a is an empirical parameter yet to be 
determined, and F is a certain functional of /(r) - the 
dimensionless ratio measuring the difference between the 
GGA exchange energy density and the exact-exchange 
energy density at a given point r. 



/(r) 



gGGA(j.)_gCxact 



(r) 



fGGA/ 



r) 



(47) 



One may note that by setting a = in Eq. (|46|) the usual 
RPA+ correction term is recovered. A simple choice 
of the functional form F{f) = f turns out to be good 
enough for fitting atomization energies, but the correct 
dissociation limit of H2 given by the standard RPA is 
destroyed. To overcome this problem, Ruzsinszky et al. 
chose a more complex form of F , 

F{f) = /[I - 7.2/2] [^ ^ i4.4/2]exp(-7.2/2), (48) 



FIG. 3: Goldstone diagrams for SOSEX contribution. The 
rules to evaluate Goldstone diagrams can be found in 
Ref. KM- 



which ensures the correct dissociation limit, while yield- 
ing significantly improved atomization energies for a — 9. 
Up to now the correction scheme of Eq. (|46p has not 
been widely benchmarked except for a small test set of 
10 molecules where the atomization energy has been im- 
proved by a factor of two [33] . 



B. Screened second-order exchange (SOSEX) 



The SOSEX correction [18|, |97|, |125[ is an important 
route to go beyond standard RPA. This concept can be 
most conveniently understood within the context of the 
ring-CCD formulation of RPA as discussed in Sec. IIIDI 
If in Eq. ((43|) . the anti-symmetrized Coulomb integrals 
Bia.jb = {ij\cib) — {ij\ba) are inserted instead of the un- 
symmetrized Coulomb integrals, the RPA-I-SOSEX cor- 
relation energy expression is obtained 



E. 



RPA+SOSEX 



1 



2 / ^ -'-ia.jb ^lajb ■ 
ij,ab 



(49) 



This approach, which was first used by Freeman l97i . 
and recently examined by Griineis et al. for solids [l25| 
and Paier [181 ] for molecular properties, has received in- 
creasing attention in the RPA community. In contrast to 
RPA 4-, this scheme has the attractive feature that it im- 
proves both total energies and energy differences simul- 
taneously. Although originally conceived in the CC con- 
text, SOSEX has a clear representation in terms of Gold- 
stone diagrams, as shown in Fig. [3] (see also Ref. [12S] ). 
which can be compared to the Goldstone diagrams for 
RPA in Fig.[2]Jc). From Fig.[3l it is clear that the leading 
term in SOSEX corresponds to the second-order exchange 
term of MP 2. In analogy, the leading term in RPA cor- 
responds to the second-order direct term of MP2. Physi- 
cally the second-order exchange diagram describes a (vir- 
tual) process in which two particle-hole pairs are created 
spontaneously at a given time. The two particles (or 
equivalently the two holes) then exchange their positions, 
and these two (already exchanged) particle-hole pairs an- 
nihilate themselves simultaneously at a later time. In 
SOSEX, similar to RPA, a sequence of higher-order dia- 
grams are summed up to infinity. In these higher-order 
diagrams, after the initial creation and exchange process, 
one particle-hole pair is scattered into new positions re- 
peatedly following the same process as in RPA, until it 
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annihilates simultaneously with the other pair at the end 
of the process. 

SOSEX is one-eleetron self-eorrelation free and ame- 
liorates the short-range over-eorrelation problem of RPA 
to a large exten t, leading to signifieantly better total en- 
ergies [93, I125J I . More importantly, the RPA underesti- 
mation of atomization energies is substantially redueed. 
However, the dissociation of covalent diatomic molecules, 
which is well described in RPA, worsens considerably as 
demonstrated in Ref. [18| and to be shown in Fig. [6l It 



was argued that the self-correlation error present in RPA 
mimics static correlation, which become s dom inant in the 
dissociation limit of covalent molecules III6I . 



Single excitation correction and its combination 
with SOSEX 



In most practical calculations, RPA and SOSEX corre- 
lation energies are evaluated using inp ut orb itals from a 
preceding KS or generalized KS (gKS) jl62l ] calculation. 
In this way both RPA and SOSEX can be interpreted as 
infinite- order summations of selected types of diagrams 
within the MBPT framework introduced in Sec. IlICi as 
is evident from Figs[2llc) and ([3|). This viewpoint is help- 
ful for identifying contributions missed in RPA through 
the aid of diagrammatic techniques. An an example, the 
second-order energy in RSPT in Eq. p8)) have contribu- 
tions from single excitations (SE) and double excitations. 
The latter gives rise to the familiar MP2 correlation en- 
ergy, which is included in the RPA-I-SOSEX scheme as 
the leading term. The remaining SE term is given by 
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where ii^^ is the self-consistent HP single-particle poten- 
tial, wmf is the mean-field potential associated with the 



reference Hamiltonian, and / 



-VV2 



the single-particle HP Hamiltonian (also known as the 
Fock operator in the quantum chemistry literature). A 
detailed derivation of Eq. ([50|) using second-quantization 
can be found in the supplemental material of Ref. [34| . 
The equivalence of Eqs. ([5T|l and (jSOj) can be readily 
confirmed by observing the relation between / and the 
single- particle reference Hamiltonian h^^: f = h^^ + 
^HF _ -MF^ g^^^ ^jjg f^p^ {iJ^\h^^\^Pa) = 0. Obviously 
for a HP reference where v^^ = -0^^, Eg. (fSTI) becomes 
zero, a fact known as Brillouin theorem |l45J |. Therefore, 
as mentioned in Section III C[ this term is not present in 
MP2 theory which is based on the HP reference. We 
note that a similar SE term also appears in 2nd-order 
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FIG. 4: Goldstone diagrams for renormalized single excitation 
contributions. Dashed lines ending with a cross denote the 
matrix element Avpq = {ipp\v^^ — v^'^^\'ipq). 



Gorling-Lev y pe rturbation theory (GL2) Il63l Il64l |. ab 
initio DFT [l6|, as well as in CC theory [l58||~IIow- 
ever, the SE terms in different theoretical frameworks 
differ quantitatively. For instance, in GL2 v^^ should 
be the exact-exchange OEP potential instead of the ref- 
erence mean-field potential. 

In Ref. [34| we have shown that adding the SE term of 
Eq. (|?T|) to RPA significantly improves the accuracy of 
vdW-bonded molecules, which the standard RPA scheme 
generally underbinds. This improvement carries over to 
atomization energies of covalent molecules and insulat- 
ing solids as shown in Ref. J126J . It was also observed 
in Ref. j34| that a similar improvement can be achieved 
by replacing the non-self-consistent HP part of the RPA 
total energy by its self-consistent counterpart. It ap- 
pears that, by iterating the exchange-only part towards 
self-consistency, the SE effect can be accounted for effec- 
tively. This procedure is termed "hybrid- RPA" , and has 
been shown to be promising even for surface adsorption 
problems j32l |. 

The SE energy in Eq. (j51[) is a second-order term in 
RSPT, which suffers from the same divergence problem 
as MP2 for systems with zero (direct) gap. To overcome 
this problem, in Ref. |3J| we have proposed to sum over 
a sequence of higher-order diagrams involving only single 
excitations. This procedure can be illustrated in terms of 
Goldstone diagrams as shown in Fig. 31 This summation 
follows the spirit of RPA and we denote it renormalized 
single excitations (rSE) [34| . The SE contribution to the 
2nd-order correlation energy in Eq. ([ST]) . represented by 
the first diagram in Fig [H constitutes the leading term 
in the rSE series. A preliminary version of rSE, which 
neglects the "off-diagonal" terms of the higher-order SE 
diagrams (by setting i ~ j = ■ ■ ■ and a = b ~ ■ ■ ■), 
was benchmarked for a tomi zation energies and reaction 
barrier heights in Ref. jl26l |. Recently we were able to 
also include the "off-diagonal" terms, leading to a refined 
version of rSE. This rSE "upgrade" does not affect the 
energetics of str ongl y bound molecules, as those bench- 
marked in Ref. [l26[. However, the interaction energies 
of weakly bound molecules improve considerably. A more 
detailed description of the computational procedure and 
extended benc hmar ks for rSE will be reported in a forth- 
coming paper [166[. However, we note that all the rSE 
results reported in Section (Sec. |V| correspond to the 
upgraded rSE. 

Diagrammatically, RPA, SOSEX and rSE are three 
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distinct infinite series of many-body terms, in which the 
three leading terms correspond to the three terms in 2nd- 
order RSPT. Thus it is quite natural to include all three 
of them, and the resultant RPA-|-SOSEX-frSE scheme 
can be viewed a renormalization of the normal 2nd-order 
RSPT. Therefore we will refer to RPA+SOSEX+rSE as 
^^renormalized second- order perturbation theory^^ or r2PT 
in the following. 



D. Other "beyond-RPA" activities 

There have been several other attempts to go beyond 
RPA. Here we will only briefly discuss the essential con- 
cepts behind these approaches without going into details. 
The interested reader is referred to the corresponding ref- 
erences. Following the ACFD formalism, as reviewed in 
Sec. IIIBi one possible route is to improve the interact- 
ing density response function. This can be conveniently 
done by adding the exchange correlation kernel (/xc) of 
time-depend ent D FT jl67l Il68| , that is omitted in RPA . 
Fuchs et al. [l69f , as well as HeBelmann and Gorling jl24l | 
have added the exact-exchange kernel to RPA, a scheme 
termed by these authors as RPA-I-X or EXX-RPA, for 
studying the H2 dissociation problem. RPA-f X or EXX- 
RPA displays a similar dissociation behavior for H2 as 
RPA: accurate at infinite separation, but slightly repul- 
sive at intermediate bond lengths. This scheme however 
gives rise to a noticeable improvement of the total en- 
ergy |l23J |. Furche and Voorhis examined the influence 
of several different local and non-local kernels on the at- 
omization energies of small molecules and the binding 
energy curves of rare-gas dimers |14l |. They found that 
semilocal /xc kernels lead to a diverging pair density at 
small inter-particle distances, and it is necessary to go to 
non-local /xc kernels to cure this. More work along these 
lines has to be done before conclusions can be drawn and 
accurate kernels become available. 

Quantum chemistry offers another route to go beyond 
the standard RPA by including higher-order exchange ef- 
fects (often termed as "RPAx" ) . There the two-electron 
Coulomb integrals usually appear in an antisymmctrizcd 
form, whereby exchange-type contributions, which arc 
neglected in standard RPA, are included automatically 
[lOlL Il70l | . The RPA correlation energy can be expressed 
as a contraction between the ring-CCD amplitudes and 
the Coulomb integrals (see Eq. (03])), or alternatively be- 
tween the coupling-strength-av erage d density matrix and 
the Coulomb integrals (see Ref. |119J ). Different flavors of 
RPA can therefore be constructed depending on whether 
one antisymmetrizes the averaged density matrix and/or 
the Coulomb integrals (sec Angyan et al. |l2l| ). Ac- 
cording to our definitions in this article, these schemes 
are categorized as different ways to go beyond the stan- 
dard RPA, while in the quantum chemistry community 
they might be simply referred to as RPA. The SOSEX 
correction, discussed in Sec. IIVBI can also be rewrit- 
ten in terms of a coupling-strcngth-avcragcd density ma- 



trix |119{ . Another interesting scheme was proposed by 
Hel3elmann [3^ , in which RPA is corrected to be exact at 
the third order of perturbation theory. These corrections 
show promising potential for the small molecules consid- 
ered in Ref. [38J. However, more benchmarks are needed 
for a better assessment. 

Both standard RPA and RPA with the exchange-type 
corrections discussed above have been tested in a r ange - 
sepa ration framework by several authors [la, [13, [2a, IllSl - 
Il22| . As mentioned briefly in the introduction, the con- 
cept of range separation is similar to the RPA -I- proce- 
dure, in which only the long-range behavior of RPA is re- 
tained. However, instead of additional corrections, here 
RPA at the short-range is completely removed and re- 
placed by semi-local or hybrid functionals. The price to 
pay is an empirical parameter that controls the range 
separation. The gain is better accuracy in describing 
molecular binding energies [ly, [l7[ , and increased compu- 
tational efficiency. The latter is due to a reduction in the 
number of required basis functions to converge the long- 
range RPA part, which is no longer affected by the cusp 
condition. More details on range -sep arated RPA can 
be found in the original references [la [TtI [28l. [TT8l - [T22 1 . 
Compared to the diagrammatic approaches discussed be- 
fore, the range-separation framework offers an alternative 
and computationally more efficient way to handle short- 
range correlations, albeit at the price of introducing some 
empiricism into the theory. 



V. APPLICATIONS 
A. Molecules 

RPA based approaches have been extensively bench- 
marked for molecular systems, ranging from the disso- 
ciat i on b e havi or of diatomic molecules [12, [13, [l3 llTI . 
I122I . I124 Il40ll . atomizatio n en ergies of small covalent 
molecules [12, [H, [Sy, [33, Il26| . interacti on energ ies of 
weakly bonded molecular comple x l28l.[33l Ill9l - ll21 1 , and 
chemical reaction barrier heights j4Cl Il26l | . The behavior 
of RPA for breaking covalent bonds was examined in its 
early days 1^ 13 1, and to day is still a topic of immense 
interest [13, Il22l Il24 Il40l | . Atomization energies of co- 
valent molecules are somewhat disappointing, because 
standard RPA, as well as its local correction (RPA-I-) , 
is not better than semi- local DFT functionals [i2[- This 
issue was subsequently referred to as "the RPA atom- 
ization energy puzzle" [36|. A solution can be fou nd in 
the beyond-RPA schemes such as RPA+SOSEX pl.[l25t 
and RPA-hSE [H, [IH, [Ull ■ Another major apphcation 
area of RPA are weakly bonded molecules. Due to the 
seamless inclusion of the ubiquitous vdW interactions, 
RPA clearly improves over conventional DFT function- 
als, including hybrids. This feature is very important 
for systems where middle-ranged non-local electron cor- 
relations play a significant role, posing great challenges 
to empirical or semi-empirical pairwisc-based correction 
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schemes. Finally, for activation energies it turned out 
that standard RPA performs remarkably well [40, Il26l | 
and the beyond RPA correction schemes that have been 
developed so far do not improve the accuracy of the stan- 
dard RPA [l26t . 

In the following, we will discuss the performance of 
RPA and its variants using representative examples to 
illustrate the aforementioned points. 



1. Dissociation of diatomic molecules 

The dissociation of diatomic molecules is an important 
test ground for electronic structure methods. The perfor- 
mance of RPA on prototypical molecules h as been exam- 
ined in a number of studies [H, [H [H, [H [HI [lH, [l2l 
|l40[. Here we present a brief summary of the behavior 
of RPA-bascd approaches based on data produced using 
our in-housc code FHI-aims |l34l ]. The numerical details 
and benchmark studies of ou r RPA implementation have 
been presented in Ref. |l33{ . In Fig. [5] the binding en- 
ergy curves obtained with PBE, MP2, and RPA-based 
methods are plotted for four molecular dimers, including 
two covalent molecules (H2 and N2), one purely vdW- 
bonded molecule (Ar2), and one with mi xed char acter 
(Be2). Dunning's Gaussian cc-pV6Z basis |l7ll . ll72| was 
used for II2 and N2, aug-cc-pV6Z for Ar2, and cc-pV5Z 
for Be2. Currently, no larger basis seems to be avail- 
able for Be2, but this will not affect the discussion here. 
Basis-set superposition errors (BSSE) are corrected us- 
ing the Boys-Bcrnardi counterpoise procedure |173l | . Also 
plotted in Fig. \5\ are accurate theoretical reference data 
for H2, Ar2, and Be2 coming respectively fro m th e full 
CI approach |l74{ . the Tang-Toennies model jl75l |. and 
the extended germinal model jl76l | . To visualize the cor- 
responding asymptotic behavior more clearly, the large 
bond distance regime of all curves is shown in Fig. [51 

RPA and RPA -I- dissociate the covalent molecules cor- 
rectly to their atomic limit at large separations, albeit 
from above after going through a positive "bump" at in- 
termediate bond distances. The fact that spin-restricted 
RPA calculations yield the correct H2 dissociation limit is 
quite remarkable, given the fact that most spin-restricted 
single-reference methods, including local and semi-local 
DFT, Hartree-Fock, as well as the coupled cluster meth- 
ods, yield an dissociation limit that is often too high in 
energy, as illustrated in Fig. [5l for PBE. MP2 fails more 
drastically, yielding diverging results in the dissociation 
limit for II2 and N2. The RPA -I- binding curves fol- 
low the RPA ones closely, with only minor differences. 
The rSE corrections are also quite small in this case, 
shifting the RPA curves towards larger binding ener- 
gies, with the consequence that the binding energy dips 
slightly below zero in the dissociation limit (see the N2 
example in Fig. [6]). This shift however leads to bet- 
ter molecular binding energies around the equilibrium 
where RPA systematically undcrbinds. The SOSEX cor- 
rection, on the other hand, leads to dramatic changes. 



Although "bump" free, RPA-I-SOSEX yields dissociation 
limits that are much too large, even larger than PBE. 
This effect carries over to r2PT at large bond distances 
where rSE does not reduce the SOSEX overestimation. 



For the purely dispersion-bonded dimer Ar2, all RPA- 
based approaches, as well as MP 2, yield the correct 
Cq/R^ asymptotic behavior, whereas the semi-local PBE 
functional gives a too fast exponential decay. Quantita- 
tively, the Cg dispersion coefficient is underesti mat ed by 
-- 9% within RPA (based on a PBE reference) [l77t . and 
SOSEX or rSE will not change this. In contrast, MP2 
overestimates the Cg value by ^ 18% [177( . Around the 
equilibrium point, RPA and RPA -I- underbind Ar2 sig- 
nificantly. The rSE correction improves the results con- 
siderably, bringing the binding energy curve into close 
agreement with the Tang-Tocnnics reference curve. The 
SOSEX correction, on the other hand, does very little in 
this case. As a consequence, r2PT resembles RPA+rSE 
closely in striking contrast to the covalent molecules. 

Be2 represents a more complex situation, in which both 
static correlation and long-range vdW interactions play 
an important role. In the intermediate regime, the RPA 
and RPA+ binding energy curves display a positive bump 
which is much more pronounced than for purely covalent 
molecules. At very large bonding distances, the curves 
cross the energy-zero line and eventually approach the 
atomic limit from below. The rSE correction moves the 
binding energy curve significantly further down, giving 
binding energies in good agreement with the reference 
values, whereas in the intermediate region a small pos- 
itive bump remains. The SOSEX correction exhibits a 
complex behavior. While reducing the bump it concomi- 
tantly weakens binding around the equilibrium distance. 
Combining the corrections from SOSEX and rSE, r2PT 
does well at intermediate and large bonding distances, 
but the binding energy at equilibrium is still noticeably 
too small. Regarding MP 2, it is impressive to observe 
that this approach yields a binding energy curve that 
is in almost perfect agreement with the reference in the 
asymptotic region, although a substantial underbinding 
can be seen around the equilibrium region. 

Summarizing this part, RPA with and without correc- 
tions shows potential, but at this point, none of the RPA- 
based approaches discussed above can produce quantita- 
tively accurate binding energy curves for all bonding sit- 
uations. It is possible, but we consider it unlikely, that 
iterating RPA to self-consistency will change this result. 
Apart from applications to neutral molecules, RPA and 
RPA-I-SOSEX studies have been carried out for the dis- 
sociation of c harged molecules. RPA fails drastically in 
this case 115l | , giving too low a total energy in the disso- 
ciation limit. Adding SOSEX to RPA fixes this problem, 
although the correction now overshoot s (with the excep- 
tion of H+) (see Ref. [H, S^, [HI [Hll for more details). 
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2. Atomization energies: the G2-I set 

One important molecular property for thermochem- 
istry is the atomization energy, given by E"^°^ — ^^ Ef^ 
where E™°^ is the ground-state energy of a molecule and 
£^f* that of the i-th isolated atom. According to this 
definition, the negative of the atomization energy gives 
the energy cost to break the molecule into its individ- 
ual atoms. Here we examine the accuracy of RPA-based 
approaches for atomization energies of small molecules. 



The RPA results for a set of 10 small organic molecules 
were reported in Furche's seminal work [12| where the un- 
derestimation of RPA for atomization energies was first 
observed. This benchmark set is included in their re- 
cent review [40|. A widely accepted represe ntati ve set 
for small organic molecules is the G2-I set |l43j . that 
contains 55 covalent molecules and will be used as an 
illustrative example here. The RPA- type results for the 
G2-I set have recently been reported in the work by Paier 
et al. [l8l.[T26t. 
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FIG. 7: Mean absolute percentage error (MAPE) for the 
G2-I atomization energies obtained with four RPA-based ap- 
proaches in addition to PBE, PBEO and MP2. 



In Fig [7] we present in a bar graph the mean absolute 
percentage error (MAPE) for the G2-I atomization ener- 
gies obtained by four RPA-based approaches in addition 
to GGA-PBE, the hybrid density-functional PBEO, and 
MP2. The actual values for the mean error (ME), mean 
absolute error (MAE), MAPE, and maximum absolute 
percentage error (MaxAPE) are listed in tabic |ll] in Ap- 
pendi x [B] The calculations were performed u sine FHI - 
aims [iMlill with Dunning's cc-pV6Z basis [l7li . [T72i . 
Reference data are taken from Ref. [l78j and corrected 



for zero-point energies. Fig [7] and Table |IT] illustrate 
that among the three traditional approaches, the hy- 
brid functional PBEO performs best, with a ME close 
to the "chemical accuracy" (1 kcal/mol = 43.4 meV). 
MP2 comes second, and PBE yields the largest error and 
shows a general trend towards overbinding. Concerning 
RPA-based approaches, standard RPA leads to ME and 
MAE that are even larger than the corresponding PBE 
values, with a clear trend of underbinding. RPA -I- docs 



not improve the atomization energies 112 1. 
in line with previous observations |12l . llSl 



All this is 



IMIill. As 



shown in the previous section, the underbinding of RPA 
for small molecules can be alleviated by adding the SO- 
SEX or rSE correction. And the combination of the two 
(i.e. r2PT) further brings the MAE down to 3.3 kcal/mol, 
comparable to the corresponding PBEO value. Whether 
the mechanism of this improvement can be interpreted in 
terms of the "multi-center non-locality of the correlation 
hole" as invoked by Ruzsinszky et al. |3a, |33] or not, is 
not yet clear at the moment. 



3. vdW interactions: S22 set 

As discussed above, one prominent feature of RPA is 
that it captures vdW interactions that arc of paramount 



importance for non-covalently bonded systems. Bench- 
marking RPA-based methods for vdW bonded systems 
is an active research field [H, [11 [13, HI [H, [H, Jsl [13, 
Tool. [Til, [in. Here we choose the S22 test set Ell as 



the illustrating example to demonstrate the performance 
of RPA for non-covalent interactions. This test set con- 
tains 22 weakly bound molecular complex of different size 
and bonding type (7 of hydrogen bonding, 8 of dispersion 
bonding, and 7 of mixed nature). Since its inception this 
test set has been widely adopted as the benchmark or 
training dataset for computational sch emes that aim at 
dealing with non-covalent interactions 
ing RPA-based approaches [2^ [Sj, [3J 
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includ- 
Il20f . The con- 
sensus emerging from these studies is as expected: RPA 
improves the binding energies considerably over semilocal 
functionals. 

Quantitatively the MAEs given by standard RPA re- 
ported for S22 by different groups show an unexpected 
spread. Specifically Zhu et al. ^ report an MAE of 2.79 
kcal/mol, Ren et al. 39 meV or 0.90 kcal/mol [SJI, and 
Eshuis and Furche 0.41 kcal/mo l [3al . The latter authors 
investigated this issue in detail jl86{ and concluded that 
the discrepancy is due to basis set incompleteness and 
BSSE. Using Dunning's correlation consistent basis sets 
plus diffuse functions and extrapolating to the complete 
basis set (CBS) limit Eshuis and Furche obtained a MAE 
of 0.79 kcal/mol with 0.02 kcal/mol uncertainty jl86| . 
These authors confirmed our observation that standard 
RPA generally underbinds weakly bound molecules. The 
basis set we have used, NAO tier 4 plus diffuse func- 
tion s fro m aug-cc-pV5Z (denoted as "fier 4 -I- a5Z-d") 
[33,[l3l, yields RPA results very close to the CBS limit. 
The results reported in Ref. [34| are, in our opinion, the 
most reliable RPA results (based on the PBE reference) 
for S22 so far. 

In Fig[8]thc relative errors (in percentage) of five RPA- 
based schemes are plotted for the molecules of the S22 
set. Results for PBE, PBEO, and MP2 are also pre- 
sented for comparison. For MP2 and RPA-based meth- 
ods, the relative errors (in percentage) for the 22 indi- 
vidual molecules are further demonstrated in Fig [9l The 
reference data were obtained using CCSD(T) and prop- 
erly extrapolated to the CBS limit by Takatani et al. 
[l87| . MP2 and RPA results are taken from Ref. ^. 
The RPA-hrSE, RPA-I-SOSEX, and r2PT results are pre- 
sented for the first time. Further details for these calcu- 
lations and an in-dep th d iscussion will be presented in 
a forthcoming paper |l66t . Figure [8] shows that PBE 
and PBEO fail drastically in this case, because these 
two functionals do not capture vdW interactions by con- 
struction, whereas all other methods show significant im- 
provement. Figure [9] further reveals that MP2 describes 
the hydrogen-bonded systems very accurately, but vastly 
overestimates the strength of dispersion interactions, par- 
ticularly for the TT-TT stacking systems. Compared to 
MP2, RPA provides a more balanced description of all 
bonding types, but shows a general trend to undcrbind. 
It has been shown that this underbinding is significantly 
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The detailed errors (ME, MAE, MAPE, and MaxAPE) 
for S22 are presented in table IIIll in Appendix [Bl Among 
the approaches we have investigated, RPA+rSE gives the 
smallest MAE for S22, while r2PT gives the smallest 
MAPE. Due to space restrictions it is not possible to in- 
clude the multitude of computational schemes that have 
emerged in recent ye ars for dealing with non-covalent 



interactions (I79l4l85 



Compared to these approaches, 
the RPA-based approaches presented here are completely 
parameter-free and systematic in the sense that they have 
a clear diagrammatic representation. Thus RPA-based 
approaches are expected to have a more general appli- 
cability, and may well serve as the reference for bench- 
marking other approaches for systems where CCSD(T) 
calculations are not feasible. 



FIG. 8: The MAPEs for the S22 test set obtained with RPA- 
based approaches in addition to PBE, PBEO, and MP2. The 
'Hier 4 + a5Z-d" basis set was used in the calculations. 
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reduced by adding SE corrections [3J]. The renormal- 
ized SE correction presented here gives rise to a more 
systematic correction to RPA, as can be seen from Fig[9l 
SOSEX shows a similar correction pattern as rSE for hy- 
drogen bonding and mixed interactions, but has little 
effect on the dispersion interaction. The r2PT scheme, 
that combines SOSEX and rSE corrections, overshoots 
for hydrogen bonding, but on average improves the de- 
scription of the other two bonding types. Figure [9] also 
reveals that tt-tt stacking configurations, as exemplified 
by the benzene dimer in the slip parallel geometry (^ 
11), represent the most challenging case for RPA-based 
methods. The relative error of RPA for this case is the 
largest. rSE provides little improvement, whereas SO- 
SEX worsens slightly. More work is needed to understand 
the origin of this failure. 



4- Reaction barrier heights 

One stringent test for an electronic structure method 
is its ability to predict chemical reaction barrier heights, 
i.e., the energy difference between the reactants and their 
transition state. This is a central quantity that dic- 
tates chemical kinetics. Semi-local density approxima- 
tions typically underestimate barrier heights [1^ Il88l |. 
RPA has already been benchmark ed fo r barrier heights 
in two independent studies [40, Il26| . Both studies 
used the test sets of 38 hydrogen-transfer barrier heights 
(HTBH38) and 38 non- hydrogen-tra nsfer barr ier heights 
(NHTBH38) designed by Zhao et al. [l8i,[l93| (together 



coined as BH76 in Ref. 1911). HTBH38 contains the for- 



ward and inverse barrier heights of 19 hydrogen transfer 
reactions, whereas NHTBH38 contains 19 reactions in- 
volving heavy atom transfers, nucleophilic substitutions, 
association, and unimolecular processes. The reference 
data were obtained using the "Weizmann-1" theory |l92l | 
- a procedure to extrapolate the CCS p(T ) results - or 
by o ther "best theoretical estimates" [l90l |. Paier et al. 
[12F 



presented results for standard RPA and "beyond 
RPA" approaches based on the PBE reference, where 
a two-point cc-pVTZ^-cc-pVQZ basis-set extrapolation 
strategy is used. In the work of Eshuis et al. [40[, stan- 
dard RPA results based on both PBE and TPSS [l93t 
refer ences were presented, where the Def2-QZVP basis 
used. 



19JI was used. The RPAOPBE resuhs for BH76 re- 
ported by both groups are very close, with an ME/MAE 
of -1.35/2.30 kcal/mol from the former and -1.65/3.10 
kcal/mol from the latter. 

The performance of RPA-based approaches, as well as 
PBE, PBEO, and MP2 for HTBH38/NHTBH38 test sets 
is demonstrated by the MAE bar graph in Fig. [TOl The 
calculations were done using FHI-aims with the cc-pV6Z 
basis set. The ME, MAE, and the maximal absolute 
error (MaxAE) are further presented in table IIVI in Ap- 
pendix [Bl In this case we do not present the relative 
errors, which turn out to be very sensitive to the compu- 
tational parameters due to some small barrier heights in 
the test set, and hence cannot be used as a reliable mca- 
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FIG. 10: The MAEs for the HTBH38 (full bars) and 
NHTBH38 (hatched bars) test sets obtained with RPA-based 
approaches in addition to PBE, PBEO, and MP2. The cc- 
pV6Z basis set was used in the calculations. 



sure of the performance of the approaches e xami ned here. 
Compared to the results reported in Ref. |l26l |. besides 
a different basis set (cc-pV6Z instead of cc-pVTZ— >-cc- 
pVQZ extrapolation), we also used the refined rSE cor- 
rection (see discussion in Section llV Cp for the RPA+rSE 
and r2PT results in Tab. IIVI which gives slightly better 
results in this case. 

On average, PBE underestimates the reaction barrier 
heights substantially, a feature that is well-known for 
GGA functionals. The hybrid PBEO functional reduces 
both the ME and MAE by more than a factor of two. 
However, the remaining error is still sizable. Standard 
RPA performs magnificently and shows a significant im- 
provement over PBEO. The performance of RPA4- is 
again v ery s imilar to standard RPA. As already noted 
in Ref. ^2^, both the rSE and the SOSEX corrections 
deteriorate the performance of RPA. This is somewhat 
disappointing, and highlights the challenge for designing 
simple, generally more accurate corrections to RPA. For- 
tunately, the errors of rSE and SOSEX are now in the 
opposite direction, and largely canceled out when com- 
bining the two schemes. Indeed, the AEs and MAEs of 
r2PT are not far from their RPA counterparts, although 
the individual errors are more scattered in r2PT as man- 
ifested by the larger MaxAE. 



B. Crystalline solids 

Crystalline solids are an important domain for RPA 
based approaches, in particular because the quantum 
chemical hierarchy of benchmark approaches cannot eas- 
ily be transfered to periodic systems. Over the years RPA 
calculatio ns have been p erfor med for a variety of sy stem s 
such as Si [li3, HIl [ill , Na [u^, h-BN [l|, NaCl [TTTI |. 



TABLE I: ME, MAE, MAPE (%), and MaxAPE (%) for the 
atomization energies (in eV/atom), lattice constants (in A), 
and bulk moduli (in G pa) of 24 crystalline solids. Results are 
taken fro m R ef. [1121 ]. The experimental atomization ener- 
gies Ref. [195( are corrected for te mpe rature effect (based on 
thermochemical correction data) [l96l | and zero-point vibra- 
tional energy. The experimental lattice constants have been 
corrected for anharmonic expansion effects. 

Atomization energies 

ME (eV) MAE (cV) MAPE (%) MaxAPE (%) 



LDA 

PBE 
RPA 
RPA-H 



-0.74 
0.15 
0.30 
0.35 



0.74 
0.17 
0.30 
0.35 



18.0 
4.5 
7.3 

8.7 



32.7 
15.4 
13.5 
15.0 



Lattice constants 



ME (A) MAE (A) MAPE (%) MaxAPE (%) 





-0.045 


0.045 


1.0 


3.7 




0.070 


0.072 


1.4 


2.7 




0.016 


0.019 


0.4 


0.9 


+ 


0.029 


0.030 


0.6 


1.1 



Bulk Moduli 



ME (GPa) MAE (GPa) MAPE (%) MaxAPE (%) 





9 


11 


9.6 


31.0 




-11 


11 


10.7 


23.7 




-1 


4 


3.5 


10.0 


+ 


-3 


5 


3.8 


11.4 



rare-gas solids [2l| , graphite |22|, |29| , and benzene crys- 
tals [231 . The most systematic benchmark study of RPA 
for crystalline solids was conducted by Harl, Schmika, 
22l . lll2{ . These authors reported "technically 



and Kresse 

converged" calculations using their VASP code and the 
projector augmented plane wave method for atomization 
energies, lattice constants, and bulk moduli of 24 repre- 
sentative crystals, including ionic compounds (MgO, LiF, 
NaF, LiCl, NaCl), semiconductors (C, Si, Ge, SiC, AIN, 
AlP, AlAs, GaN, GaP, GaAs, InP, InAs, InSb), and met- 
als (Na, Al, Cu, Rh, Pd, Ag). The error analysis of their 
RPA and RPA-t- results, based on a PBE reference, as 
well as the LDA and PBE results are presented in Tab.[Tl 
As is clear from Tab. HI the RPA lattice constants and 
bulk moduli are better than in LDA and PBE. The at- 
omization energies, however, are systematically underes- 
timated in RPA, and the MAE in this case is even larger 
than that of PBE. This behavior is very similar to that for 
the atomization energies in the G2 set discussed above. 
Harl et al. also observed that the error of RPA does not 
grow when going to heavier atom s, or open-shell systems 
in contrast to LDA or PBE flT^ . The RPA-f results are 
in general slightly worse than those of standard RPA. 

The performance of RPA has not been extensively 
benchmarked for vdW-bound solids. However, judging 
from the studies on h-BN [1^, ra re-g as crystals |2l| . 
benzene crystal [2j], and graphite [2^, standard RPA 
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does an excellent job regarding the equilibrium lattice 
constants and cohesive energies, whereas semi-local DFT 
fails miserably, yielding typically too weak binding and 
too large lattice constants. LDA typically gives a finite 
binding (often overbinding) for weakly bonded solids, but 
its performance varies significantly from system to sys- 
tem, and cannot be trusted in general, as the functional 
by construction does not contain the necessary physics 
to reliably describe this phenomenon. In a recent work 
^^ , it was shown that RPA also reproduces the correct 
1/d^ asymptotics between g raph ite layers as analytically 



predicted by Dobson et al. jl97| . This type of behavior 



can neither be described by LDA, GGAs, nor by hybrid 
functionals. 

For solids attempts have also been made to go beyond 
the standar d RP A. For a smaller test set of 11 insulators, 
Paier et al. |l26| showed that adding SOSEX corrections 
to RPA the MAE of atomization energies is reduced from 
0.35 eV/atom to 0.14 eV/atom. By replacing the non- 
self-consistent HF energy by its self-consistent counter- 
part, which mimicks the effect of adding single excita- 
tion corre ction s [3J|, reduces the MAE further to 0.09 
eV/atom jl26| . Thus the trend in periodic insulators is 
again in line with what has been observed for molecular 
atomization energies. The effects of SOSEX and rSE cor- 
rections for metals, and for other properties such as the 
lattice constants and bulk moduli have not been reported 
yet. 



C. Adsorption at surfaces 

The interaction of atoms and molecules with surfaces 
plays a significant role in many phenomena in surface 
science and for industrial applications. In practical cal- 
culations, the super cells needed to model the surfaces 
are large and a good electronic structure approach has 
to give a balanced description for both the solid and the 
adsorbate, as well as the interface between the two. Most 
approaches today perform well for either the solid or the 
isolated adsorbate (e.g. atoms, molecules, or clusters), 
but not for the combined system, or are computationally 
too expensive to be applied to large super cells. This is 
an area where we believe RPA will prove to be advanta- 
geous. 

The systems to which RPA has been applied include Xe 
and 3,4,9, 10-perylene-tetracarbox ylic acid dianhydride 
(PTCDA) adsorbed on Ag(lll) [25|; CO on Cu(lll) 
[23 . l26j and other noble/transiti on m etal surfaces [27| : 
benzene on Ni(lll) ^, Si(OOl) [lil|, and the graphite 
surface jll3{ : and graphene on N i(lll) jl98l Il99f . and 
Cu(lll), Co(OOOl) surfaces ^9^. In all these applica- 
tions, RPA has been very successful. 

To illustrate how RPA works for an adsorbate system, 
here we briefly describe the RPA study of CO@Cu(lll) 
following Ref. [2^. The work was motivated by the so- 
called "CO adsorption puzzle" - LDA and several GGAs 
predict the wrong adsorption site for CO adsorbed on 
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FIG. 11: Adsorption energies for CO adsorbed at the on- 
top and fee hollow sites of the Cu(lll) surface as obtained 
using LDA, AM05, PBE, PBEO, and RPA. RPA results are 
presented for both PBE and PBEO references, and they differ 
very little. 



sever al noble/transition metal surfaces at low coverage 
[200| . For instance, for the (111) surface of Cu and Pt 
DFT within local/semi-local approximations erroneously 
favor the threefold-coordinated hollow site, whereas ex- 
periments clearly show that the singly-coord inated on - 
top site is the energetically most stable site [20 ll |202| | . 
This posed a severe challenge to the first-principles mod- 
eling of inolecular adsorption problems and the question 
arose, at what level of approximation can the correct 
physics be recovered. In our study, the Cu surface was 
modeled using systematically increasing Cu clusters cut 
out of the Cu(lll) surface. Following a procedure pro- 
posed by Hu et al. J203l ] , the RPA adsorption energy was 
obtained by first converging its difference to the PBE val- 
ues with respect to cluster size, and then adding the con- 
verged difference to the periodic PBE results. The RPA 
adsorption energies for both the on-top and fee (face cen- 
tered cubic) hollow sites are presente d in Fig. II 11 together 
with the resuhs from LDA, AMOS ^0^, PBE, and the 
hybrid PBEO functional. Fig.lTTjreveals what happens in 
the CO adsorption puzzle when climbing the so-called Ja- 
cob's ladder in DFT (6[ — going from the first two rungs 
(LDA and GGAs) to the fourth (hybrid functionals), and 
finally to the fifth rung (RPA and other functionals that 
exphcitly depend on unoccupied KS states). Along the 
way the magnitude of the adsorption energies on both 
sites are reduced, but the effect is more pronounced for 
the fee hollow site. The correct energy ordering is already 
restored at the PBEO level, but the energy between the 
two sites is too small. RPA not only gives the correct 
adsorption site, but also produces a reasonable adsorp- 
tion energy difference of 0.22 eV, consistent with exper- 
iments. This result was later confirmed by the periodic 
RPA calculations of Harl and Kresse in Ref. [23, with 
only small numerical differences arising from the differ- 
ent implementations and different convergence strategy. 
The work of Rohlfing and Brcdow on Xe and PTCDA 
adsorbed at Ag(lll) surface represents the first RPA 
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study regarding surface adsorption problems, where the 
authors exphcitly demonstrated that RPA yields the ex- 
pected — C3/((i — do)^ behavior for large molecule-surface 
separations d. Schimka et al. extended the RPA bench- 
mark studies of the CO adsorption problem to more sur- 
faces [27[ . They found that RPA is the only approach so 
far that gives both good adsorption energies as well as 
surface energies. GGAs and hybrid functionals at most 
yield either good surface energies, or adsorption ener- 
gies, but never both. Goltl and Hafner investigated the 
adsorption of small alkanes in Na-exchanged chabazite 
using RPA and several other approaches. They found 
that the "hybrid RPA" scheme, as p roposed in Ref. [3J| 
and further examined in Ref. jl26l |. provides the most 
accurate description of the system compare d to the al- 
ternatives e.g. DFT-D [lol and vdW-DF ^^. More 
recently RPA was applied to the adsorption of benzene 
on the Si(OOl) surface by Kim et al |ll 4|. gr aphcne on the 
Ni(lll ) sur face by Mittendorfer et al |l98j | and by Olsen 
et al. |l99t . and additionally graphene on Cu(lll) and 
Co(OOOl) surfaces by the latter authors. In all these stud- 
ies, RPA is able to capture the delicate balance between 
covalcnt and dispersive interactions, and yields quantita- 
tively reliable results. We expect RPA to become increas- 
ingly more important in surface science with increasing 
computer power and more efficient implementations. 



VI. DISCUSSION AND OUTLOOK 

RPA is an important concept in physics and has a more 
than 50 year old history. Owing to its rapid develop- 
ment in recent years, RPA has shown great promise as a 
powerful first-principles electronic-structure method with 
significant implications for quantum chemistry, compu- 
tational physics, and materials science in the foreseeable 
future. The rise of the RPA method in electronic struc- 
ture theory, and its recent generalization to r2PT, were 
borne out by realizing that traditional DPT functionals 
(local and semi-local approximations) are encountering 
noticeable accuracy and reliability limits and that hybrid 
density functionals are not sufficient to overcome them. 
With the rapid development of computer hardwares and 
algorithms, it is not too ambitious to expect RPA-based 
approaches to become (or at least to inspire) main-stream 
electronic-structure methods in computational materials 
science and engineering in the coming decades. At this 
point it would be highly desirable if the community would 
start to build up benchmark sets for materials science 
akin to t he ones in quantum chemistry (e.g. G2 |143| or 
S22 [l44| ). These should include prototypical bulk crys- 
tals, surfaces, and surface adsorbates and would aid the 
development of RPA-based approaches. 

As an outlook, we would like to indicate several direc- 
tions for future developments of RPA-based methods. 

i. Improved accuracy: Although RPA does not suffer 
from the well-documented pathologies of LDA and 
GGAs, its quantitative accuracy is not always what 



is desired, in particular for atomization energies. 
To improve on this and to make RPA worth its 
computational effort, further corrections to RPA 
are necessary. To be useful in practice, these should 
not increase the computational cost significantly. 
The r2PT approach as presented in Sec. IIVI and 
benchmarkcd in Sec. |V]is one example of this kind. 
More generally the aim is to develop RPA-based 
computational schemes that are close in accuracy 
to CCSD(T), but come at a significantly reduced 
numerical cost. More work can and should be done 
along this direction. 

ii. Reduction of the computational cost: The major 
factor that currently prevents the widespread use 
of RPA in materials science is its high numeri- 
cal cost compared to traditional DFT methods. 
The state-of-the-art implementations still have an 
0{N'^) scaling, as discussed in Sec. IIIII To enlarge 
the domain of RPA applications, a reduction of this 
scaling behavior will be highly desirable. Ideas can 
be borrowed from 0{N) methods [l39l | developed 
in quantum-chemistry (in particular in the context 
of MP 2) or c omp ression techniques applied in the 
GW context [T38l |. 

iii. RPA forces: For a ground-state method, one cru- 
cial component that is still missing in RPA are 
atomic forces. Relaxations of atomic geometries 
that are common place in DFT and that make DFT 
such a powerful method are currently not possible 
with RPA or at least have not been demonstrated 
yet. An efficient realization of RPA forces would 
therefore extend its field of application to many 
more interesting and important materials science 
problems. 

iv. Self- consistency: Practical RPA calculations are 
predominantly done in a post-processing manner, 
in which single-particle orbitals from KS or gen- 
eralized KS calculations are taken as input for a 
one-shot RPA calculation. This introduces undc- 
sired uncertainties, although the starting-point de- 
pendence is often not very pronounced, if one re- 
stricts the input to KS orbitals. A self-consistent 
RPA approach can be defined within KS-DFT 
via the optimized effective potential method j95| . 
and ha s already been applied in a few instances 
SH, El [Hi, [ila, [HI]. However, in its current re- 
alizations self-consistent RPA is numerically very 
challenging, and a more practical, robust and nu- 
merically more efficient procedure will be of great 
interest. 

With all these developments, we expect RPA and its 
generalizations will play an increasingly important role 
in computational materials science in the near future. 
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Appendix A: RI-RPA implementation in FHI-aims 



where 



In this section we will briefly desc ribe how RPA 
is implemented in the FHI-aims code |134| using the 
resolution-of-ide ntity (RI) technique. More details can 
be found in R ef. [l33l . For a different formulation of RI- 
RPA see Ref. 30, So] • We start with the expression for 
the RPA correlation energy in Eq. (|23|) . which can be 
formally expanded in a Taylor series, 

1 r°° °° 1 



E. 



RPA 



duj > — Tr 

n=2 



x"(i^)v) 



(Al) 



Applying RI to RPA in this context means to represent 
both x°(*w) and v in an appropriate auxiliary basis set. 
Eq. (jAip can then be cast into a series of matrix op- 
erations. To achieve this we perform the following RI 
expansion 



V'*(r)V',(r) 






C^.P,iv). 



(A2) 



where P^{v) are auxiliary basis functions, C^, are the ex- 
pansion coefficients, and A'aux is the size of the auxiliary 
basis set. Here C serves as the transformation matrix 
that reduces the rank of all matrices from iVocc * -^vir 
to Aaux, with A'occ, ^vir and A^aux being the number 
of occupied single-particle orbitals, unoccupied (virtual) 
single-particle orbitals, and auxiliary basis functions, re- 
spectively. The determination of the C coefficients is 
not unique, but depends on the underlying metric. In 
quantum chemistry the "Coulomb metric" is the stan- 
dard choice where the C coefficients are determined by 
minimizing the Coulomb repulsion between the resi duals 
of the expansion in Eq. (|A2p (for details see Ref. [l33J | 
and references therein). In this so-called "RI-V" approx- 
imation, the C coefficients are given by 



where 



and 
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0.(r)0,(r)P,(r') 



drdr' . 



V, 
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^^^m^drdr'. 
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(A4) 



(A5) 



In practice, sufficiently accurate auxiliary basis set can 
be constructed such that N^u^ <^ Nocc * ^vir, thus reduc- 
ing the computational effort considerably. A practically 
accurate and efficient way of constructing auxiliary basis 
set {P^{r)} and their associated {C^A for atom-centered 
basis functions of general shape has been presented in 
Ref. [133]. 

Combining Eq. ^ with (|A2|) yields 
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Introducing the Coulomb matrix 



V^,^ 



JJdrr'P^{r)v{r,r')P,{r'). 
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we obtain the first term in (jAip 
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= - ^ / duj Y^ x°,.(Jw)K„x°/3(iw)F^, 






(A9) 



This term corresponds to the 2nd-order direct correla- 
tion energy also found in MP2. Similar equations hold 
for the higher order terms in (|Aip . This suggests that 
the trace operation in Eq. ([23|) can be re-interpreted as a 
summation over auxiliary basis function indices, namely, 
Tr[AS] = T.^,u^^^'^B'^^^■- provided that x°(r,r',iw) and 
v{r, r') are represented in terms of a suitable set of aux- 
iliary basis functions. Equations (P5)) . (jA7p . and (|A8[) 
constitute a practical scheme for RI-RPA. In this imple- 
mentation, the most expensive step is Eq. ()A7p for the 
construction of the independent response function, scal- 
ing as N'^^^N occNvir- We note that th e sam e is true for 
standard plane-wave implementations |112| where iVaux 
corresponds to the number of plane waves used to expand 
the response function. In that sense RTbased local-basis 
function implementations are very similar to plane- wave- 
based or LAPW-based implementations, where the plane 
waves themselves or the mixed product basis serve as the 
auxiliary basis set. 



Appendix B: Error statistics for G2-I, S22, and 
NHBH38/NHTBH38 test sets 

Tables [Hi IIIH and IIVI present a more detailed error 
analysis for the G2-I, S22, and NHBH38/NHTBH38 test 
sets. Given are the mean error (ME), the mean ab- 
solute error (MAE), the mean absolute percentage er- 
ror (MAPE), the maximum absolute percentage error 
(MaxAPE). and the maximum absolute error (MaxAE). 
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TABLE II: ME (in eV), MAE (in eV), MAPE (%), 
MaxAPE(%) for atomization energies of the G2-I set obtained 
with four RPA-based approaches in addition to PBE, PBEO 
and MP2. A negative ME indicates overbinding (on average) 
and a positive ME underbinding. The cc-pV6Z basis set was 
used. 





ME 


MAE 


MAPE 


MaxAPE 


PBE 


-0.28 


0.36 


5.9 


38.5 


PBEO 


0.07 


0.13 


2.6 


20.9 


MP2 


-0.08 


0.28 


4.7 


24.6 


RPA 


0.46 


0.46 


6.1 


24.2 


RPA+ 


0.48 


0.48 


6.3 


24.9 


RPA+SOSEX 


0.22 


0.25 


4.2 


36.7 


RPA+rSE 


0.30 


0.31 


4.0 


22.5 


r2PT 


0.07 


0.14 


2.6 


24.9 



TABLE III: ME (in meV), MAE (in meV), MAPE (%), and 
MaxAPE(%) for the S22 test set [iti] obtained with five RPA- 
based approaches in addition to PBE, PBEO, and MP2 ob- 
tained with FHI-aims. The basis set ^Hier 4 -I- a5Z-d" [l33l ] 
was used in all calculations. 





ME 


MAE 


MAPE MaxAPE 


PBE 


116.2 


116.2 


57.8 


170.3 


PBEO 


105.7 


106.5 


55.2 


169.1 


MP2 


-26.5 


37.1 


18.7 


85.1 


RPA 


37.8 


37.8 


16.1 


28.7 


RPA+ 


51.2 


51.2 


21.9 


39.4 


RPA-hrSE 


14.1 


14.8 


7.7 


24.8 


RPA-hSOSEX 


15.4 


18.0 


10.5 


34.6 


r2PT 


-8.4 


21.0 


7.1 


30.7 



TABLE IV: ME, MAE, and MaxAE (in eV) for the HTBH38 
[l89l | and NHTBH38 ^^ test sets obtained with four RPA- 
based approaches in addition to PBE, PBEO, and MP2, as 
obtained using FHI-aims. The cc-pV6Z basis set was used in 
all calculations. Negative ME indicates an underestimation 
of the barrier height on average. 





HTBH38 


NHTBH38 




ME MAE MaxAE 


ME MAE MaxAE 


PBE 


-0.399 0.402 


0.863 


-0.365 0.369 


1.320 


PBEO 


-0.178 0.190 


0.314 


-0.134 0.155 


0.609 


MP2 


0.131 0.169 


0.860 


0.215 0.226 


1.182 


RPA 


0.000 0.066 


0.267 


-0.065 0.081 


0.170 


RPA+ 


0.005 0.069 


0.294 


-0.068 0.084 


0.168 


RPA-frSE 


-0.170 0.187 


0.809 


-0.251 0.252 


0.552 


RPA+SOSEX 


0.243 0.244 


0.885 


0.185 0.188 


0.781 


r2PT 


0.072 0.084 


0.453 


-0.001 0.129 


0.432 
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